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ABSTRACT 

We present timing observations of four millisecond pulsars discovered in the Parkes 
20-cm multibeam pulsar survey of the Galactic plane. PSRs J1552—4937 and 
J1843—1448 are isolated objects with spin periods of 6.28 and 5.47 ms respectively. 
PSR J1727—2946 is in a 40-day binary orbit and has a spin period of 27 ms. The 
4.43-ms pulsar J1813—2621 is in a circular 8.16-day binary orbit around a low-mass 
companion star with a minimum companion mass of 0.2 Mq. Combining these re¬ 
sults with detections from five other Parkes multibeam surveys, gives a well-defined 
sample of 56 pulsars with spin periods below 20 ms. We develop a likelihood analysis 
to constrain the functional form which best describes the underlying distribution of 
spin periods for millisecond pulsars. The best results were obtained with a log-normal 
distribution. A gamma distribution is less favoured, but still compatible with the obser¬ 
vations. Uniform, power-law and Gaussian distributions are found to be inconsistent 
with the data. Galactic millisecond pulsars being found by current surveys appear 
to be in agreement with a log-normal distribution which allows for the existence of 
pulsars with periods below 1.5 ms. 
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1 INTRODUCTION 

Millisecond radio pulsars (MSPs) are fascinating objects to 
study. Their phenomenal rotational stability allows them 
to be used for a wide variety of fundamental physics 
experiments including as a Galactic-scale observatory to 
search for low-frequency gravitational waves (see, e.g.. 
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iHobbs et al.ll2010l). Ever since the discovery of the first MSP 
llBacker et al.l Il982l ) it has been clear that the difficulties 
in detection imply that the Galactic population of MSPs 
is substantial. Early studies showed that the population of 
MSPs is comparable to that of normal pulsars ( see,.e.g., 
iKulkarni fc Naravanlll98^ : Ijohnston fc Baii^IlQQlI) . 

The continued improvement of data acquisition systems 
over the past twenty years has led to a dramatic increase in 
survey sensitivity to MSPs. The number of known MSPs in 
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Figure 1. Integrated pulse profiles each showing 360 degrees of rotational phase at 20 cm wavelength for the four MSPs described in 
this paper. Data were obtained with the Parkes digital filterbank systems. 
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Figure 2. Integrated polarization profiles for the two MSPs for which significant RMs were detected. The top panels of each plot show 
the polarization position angle as a function of pulse phase. The bottom panels show total intensity (bold line), linear polarization (solid 
line) and circular polarization (dotted line). 


the Galactic disk (i.e ., those not associated with globular 
clusters) is now 23 QtL Because of the great su ccess of blind 
surve ys of the Galactic field (for a review, see IStovall et al ' 


2013), and targeted searches of Fermi sources (IRav et al 


20121 ). for the first time in a decade, Galactic MSPs out¬ 


number their counterparts in globular clusters. 

The Parkes multibeam pulsar survey (PMPS) of the 
Galactic plane is the most successful large-scale search for 
pulsars so far undertaken. Six previous papers in this se¬ 
ries have presented timing parameters for 742 newly dis- 


survey results (IManchester et al.l 200ll: Morris et al. 

2 OO 2 I; 

Kramer et al. 

I 2 OO 3 I: Hobbs et al.ll2004l; IFaulkner et al.l 

20o4 

Lorimer et al. 

20061: iGrawford et al. l2013l). Over the nast 


five years, severa l re-analyse s of t he survey data have 
been carried out. iKeith et al.l (l2009l) discovered a further 
28 pulsars by applying new candi date sorting alg o rithm s 
to the data process ed earlier by iFaulkner et ahl ||2004|) . 
lEatough et al.l ll2009h applied a new interference removal 
technique to a small portion of the data and discovered 
a further four pulsars. In a fu rther reanalysis, K eane et 
al. found one fast radio b urst (IK eane et al.l 2013) and 10 
rotating radio transients dKeane et al. 201C ) in a dditio n 
to the 11 foun d orig inally by iMcLauehlin et al.l ll2006l) . 
iMickaliger et al.l ll2012r) reported the discovery of the 34.5 ms 
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binary pulsar J1725—3853 as well as four other millisecond 
pulsars. One other binary MSP, J1753—2814, has also been 
discovered as a result of this processing effort (Mickaliger 
et al. in preparation). Following e arlier discoverie s usin g 
“stack-slide” a c celera tion searches bv iFaulkner et aP ll2004h . 
lEatough et al.l (l2013f) report the discovery of 16 pulsars in 
a coherent acceleration search of th e data. Ongoing pro cess- 
ing by Einstein@Home volunteers llKnisoel et al.ll2013l) has 
resulted in the discovery of a further 23 pulsars. 


In this paper, we present timing solutions for four MSPs 
discovered in the PMPS. Preliminary discovery and confir¬ 
mation o bservations of these p ulsars were previously pub¬ 
lished by IFaulkner et ^ (|2004| ). The total number of pul¬ 
sars found in the survey so far stands at 833. Since ex- 
tensiye population studies of the normal pulsar popula¬ 
tion as reyealed by the PMPS haye already been car ried 
out llFaucher-Gigufere fc Kaspill2006l : lLorimer et al.ll2006l ). in 
this paper we focus our discussion on the spin period dis¬ 
tribution of the MSP population. The plan for this paper 
is as follows. In iJ5]we present the basic timing parameters, 
pulse widths, mean profiles and flux densities for the four 
new MSPs. In ^we compile a sample of MSPs and use it 
to carry out a likelihood analysis to constrain the underly¬ 
ing distribution of spin periods. The main conclusions are 
summarized in Sjd] 
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Table 1. Spin, astrometric and derived parameters from the timing analysis of four MSPs. All astrometric parameters are given in the 
J2000 coordinate system. The reduced values from each fit is listed as x?- Figures in parentheses represent the uncertainties in the 
least significant digit and are the nominal l-tr TEMP02 uncertainties. Distance estimates are based on the pulsar DM using the Cordes 
& Lazio (2002) NE2001 electron-density model. Pseudo-luminosities are computed by multiplying flux density by distance squared. 
Characteristic ages, magnetic fields and spin-down luminosities are based on the spin period and period derivative (see, e.g., Lorimer & 
Kramer 2005) and account for the contributions due to the Shklovskii effect and Galactic acceleration. 


Parameter 

PSR J1552-4937 

PSR J1727-2946 

PSR J1813-2621 

PSR J1843-1448 

R.A. (hh:mm:ss.s) 

15:52:13.2709(4) 

17:27:15.09493(17) 

18:13:40.59165(10) 

18:43:01.3750(3) 

Dec. (dd:mm:ss.s) 

-49:37:49.744(11) 

-29:46:36.797(17) 

-26:21:57.055(18) 

-14:48:12.61(3) 

Proper motion in R.A. (mas y“^) 

-3(3) 

0.6(9) 

-7.3(9) 

10.5(19) 

Proper motion in Dec. (mas y~^) 

-13(8) 

0(8) 

-22(16) 

12(15) 

Epoch of position 

J2000 

J2000 

J2000 

J2000 

Spin period (ms) 

6.2843113814174(12) 

27.0831832440066(12) 

4.4300116286341(18) 

5.4713308755095(6) 

First derivative of spin period 

1.900(4) X 10-20 

2.4632(3) X 10-1® 

1.2466(6) X 10-2° 

6.209(18) X 10-21 

Dispersion measure (cm“^ pc) 

114.19(8) 

60.74(3) 

112.524(9) 

114.51(7) 

Rotation measure (rad m~^) 

— 

-61(32) 

136(8) 

— 

Epoch of period (MJD) 

54033 

54723 

54058 

53934 

Data span (MJD) 

52860-55206 

52666-56781 

52696-55419 

52696-55152 

/ degrees of freedom 

1.00/132 

1.45/208 

0.90/134 

1.16/115 

Post-fit rms residual (ps) 

78 

43 

17.5 

49 

Flux density at 1.4 GHz (mjy) 

0.14 

0.25 

0.65 

0.57 

Pulse width at 50% of peak (ms) 

0.9 

1.8 

0.66 

1.0 

Distance (kpc) 

4.8 

1.4 

2.9 

2.9 

Pseudo-luminosity (mJy kpc^) 

3.2 

0.49 

2.1 

4.8 

Intrinsic period derivative 

0.7(1.1) X 10-21 

2.426(7) X 10-1® 

-0.6(1.7) X 10-21 

-0.5(1.1) X 10-21 

Characteristic age (Gy) 

> 15 

1.77 

> 5.6 

> 14 

Surface magnetic field (10® G) 

< 2.1 

25.9 

< 2.4 

< 1.9 

Spin-down luminosity (10®® ergs s“^) 

< 1.1 

0.48 

< 5.7 

< 1.5 


2 FOUR MILLISECOND PULSARS 


The pulsars were discovered using the processing schemes 
described by iFaulkner et al.l (|2004) . Following the confir- 
mation and p ositio nal refinement procedures described by 
I Morris et al.l (l2002l ). each pulsar was observed regularly 
at Parkes using initially the 512 x 0.5 MHz analogue fil- 
terbank system ([Manchester et ^ 2001 ) and subsequ ently 
the digital filterbank systems ( Manchester et ^l2013l ). For 
each pulsar, pulse times of arrival were determined from 
the individual obse rvations using standard pulsar timing 
techniques (see, e.g.. lLorimer fc Krameill 2005l) im pleme nted 
in the psrchive software package ( Hot an et al. | [2004lfl . A 
model containing the spin, astrometric and (if necessary) 
any binary parameters was fitt ed to the arrival t imes us¬ 
ing the TEMP02 timing package dHobbs et al.ll20(3^ ). Arrival 
times were referred to TT(TA I) and the DE421 planetary 
ephemeris (iFolkner et al.l 12^0081 ) was used. Timing parame- 
ters are expressed in “TCB” units native to tempo2 (see 
iHobbs et al.l l2006l . for the definition of TCB). Integrated 
pulse profiles are shown in Fig. [T] 

Timing parameters from these analyses along with var¬ 
ious derived quantities are presented in Tables [T] and O For 
PSR J1727-2947, time-of-arrival uncertainties were multi¬ 
plied by a factor ranging between 0.85-1.5 for different back¬ 
end systems to maintain a reduced value close to unity. 
Also listed in Table [T] is the post-fit root-mean-square resid¬ 
ual. The values obtained from our timing so far are relatively 
large (17-78 /rs) and indicate that these pulsars are un- 
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likely to be useful additions to MSP timing arrays. Although 
proper motions in right ascension have been measured for 
PSRs J1813—2621 and J1843—1448, we are unable to mea¬ 
sure a significant proper motion in declination because of 
the low ecliptic latitude of these pulsars. Flux densities at 
1400 MHz and pulse widths at 50% of the peak level based 
on the profiles shown in Fig. [T]are listed in Table [T] 

For two of the MSPs, J1727-2947 and J1813-2621, sig¬ 
nificant levels of polarized emission was measured and these 
are shown in the integrated pulsed profiles in Fig. (2] Rota¬ 
tion measures for both these pulsars were determined using 
the rmfit tool within psrchive with conservative estimates 
of the uncertainties. 

PSRs J1552—4937 and J1843—1448 bring the total 
number of isolated MSPs known in the Galactic disk to 
37. When compared to the sample of 172 MSPs for which 
an orbiting companion has been confirmed, the fraction of 
observed isolated MSPs currently stands at 18%. An out¬ 
standing issue in our understanding of MSP population is 
to explain this population in a self-consistent fashion. In par¬ 
ticular, an open question is whether isolated MSPs formed 
in a different way from binary MSPs. We discuss this issue 
further in § 13.61 

PSR J1727—2947 is a relatively long-period MSP [P ~ 
27 ms) in a mildly eccentric (e ~ 0.04) 40-day binary system. 
With a minimum companion mass of ~ 0.8 Mq, the system 
is most likely a member of th e so-called “interm ediate-mass 
binary pulsar” (IMBP) class (|Camilo et al.lf200ll ) with a rel¬ 
atively massive CO white dwarf companion. The parame¬ 
ters for PSR J1813—2621 imply that it is very represen¬ 
tative of the low-mass binary MSP population. Interpret¬ 
ing the orbital parameters in the standard way (see, e.g.. 
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Table 2. Measured and derived orbital parameters for 
PSRs J1727—2946 and J1813—2621 which use the “BT” binary 
model telandford &: TeukolskAl Il976h and the “ELLl” binary 
model (Lange et al. 2001) respectively. Parameters listed are the 
binary period (Ph), projected semi-major axis (asini), orbital 
eccentricity (e), first and second Laplace-Lagrange parameters 
(ei and 62 ), longitude and epoch of periastron {u and Tq) and 
epoch of ascending node (Tasc)- Figures in parentheses represent 
l-(j TEMP02 uncertainty in the least significant digits. For PSR 
J1813—2621, we also list the corresponding values of e, lj and 
To computed from the Laplace-Lagrange parameters. The Keple- 
rian mass function (An^G(asmi)^/P^, where G is Newton’s con¬ 
stant), and the minimum companion mass (calculated assuming 
a 1.4 Mq pulsar and setting, i = 90°) are listed. 


PSR 

J1727-2946 

J1813-2621 

Binary model 

BT 

ELLl 

Pb (d) 

40.30771094(3) 

8.159760702(10) 

al sin 2 (It sec) 

56.532497(5) 

5.592583(3) 

e 

0.04562943(16) 


^ n 

320.39625(20) 

289^1® 

To (MJD) 

54711.47169(3) 

54061.51®° 

ei 

- 

-2.5(10) X 10“° 

e2 

- 

9(8) X 10-'^ 

Tasc (MJD) 

- 

54054.9328319(6) 

Mass function (Mq) 

0.1194 

0.00282 

Min. comp, mass (Mq) 

0.827 

0.188 


iLorimer fc Kram^l2005l) . we infer a companion mass of at 
least 0.2 Mq, typical of a low-mass white dwarf. Optical 
studies of these companions may provide further insights 
into the nature of these two binary systems. 

For nearby MSP s, it is well known (see, e.g., 
iDamour fc Tavlol^ll991^ that two significant contributions 
to the observed period derivative are the effects of secular ac¬ 
celeration (sometimes referred to as the “Shklovskii effect”, 
IShklovs^IlQTOl i and Galactic acceleration. For a pulsar of 
period P, trausverse speed V, acceleration a, distance D 
with an intrinsic period derivative Pint, the observed period 
derivative 


Pobs — Pint + P 



( 1 ) 


where n is a unit vector along the line of sight to the pulsar 
and c is the speed of l ight. Following the discussion in §3.1 of 
I Nice &: Tavloil (Il995h to compute these effects, and assum¬ 
ing a 25% uncertainty on the di stances, computed usin g the 
NE2001 electron density model ICordes fc Laziol (l2002l l. we 
calculated or placed limits on Pint for each pulsar and list 
our results in Table 1. The resulting characteristic age and 
magnetic field strength estimates for these pulsars are also 
indicated in Table 1. As can be seen, for all but PSR J1727- 
2946, these effects account for most of the observed period 
derivative. 


3 THE SPIN PERIOD DISTRIBUTION OF 
GALACTIC MSPS 

The large sample of over 1000 normal pulsars detected in 
the various Parkes multibeam surveys has provided signif¬ 
icant advances in our understanding of the normal pul- 


Table 3. The 56 MSPs used in the population study. For each 
pulsar, we list the spin period (P), dispersion measure (DM), 
Galactic longitude (Z), Galactic latitude (b), whether this is a 
binary pulsar as well as the survey which detected the pulsar. 
The surveys considered were the PMPS, the Parkes high-latitude 
(PH) pulsar survey, the Perseus arm (PA) pulsar survey, the Deep 
Multibeam (DMB) survey, the Swinburne intermediate latitude 
(SWIL) and high latitude (SWHL) surveys. 


PSRJ 

P 

(ms) 

DM 

(pc/cc) 

1 

n 

b 

n 

Bin 

Survey 

0437-4715 

5.76 

2.6 

253.4 

-42.0 

Y 

PH 

0610-2100 

3.86 

60.7 

227.7 

-18.2 

Y 

PH 

0711-6830 

5.49 

18.4 

279.5 

-23.3 

N 

SWHL 

0721-2038 

15.54 

76.1 

234.7 

-2.9 

Y 

PA 

0900-3144 

11.11 

75.7 

256.2 

9.5 

Y 

PH 

0922-52 

9.68 

122.4 

273.8 

-1.4 

Y 

PMPS 

1022-1-1001 

16.45 

10.2 

231.8 

51.1 

Y 

PH 

1024-0719 

5.16 

6.5 

251.7 

40.5 

N 

PH 

1045-4509 

7.47 

58.2 

280.9 

12.3 

Y 

SWIL 

1125-6014 

2.63 

53.0 

292.5 

0.9 

Y 

PMPS 

1147-66 

3.72 

133.5 

296.5 

-4.0 

Y 

PMPS 

1216-6410 

3.54 

47.4 

299.1 

-1.6 

Y 

PMPS 

1435-6100 

9.35 

113.7 

315.2 

-0.6 

Y 

PMPS 

1546-59 

7.79 

168.2 

323.5 

-3.8 

Y 

PMPS 

1552-4937 

6.28 

114.6 

330.0 

3.5 

N 

PMPS 

1600-3053 

3.60 

52.3 

344.1 

16.5 

Y 

SWIL 

1603-7202 

14.84 

38.0 

316.6 

-14.5 

Y 

SWIL 

1618-39 

11.99 

117.5 

340.8 

7.9 

Y 

SWIL 

1629-6902 

6.00 

29.5 

320.4 

-13.9 

N 

SWIL 

1643-1224 

4.62 

62.4 

5.7 

21.2 

Y 

SWHL 

1652-48 

3.78 

187.8 

337.9 

-2.9 

Y 

PMPS 

1708-3506 

4.50 

146.8 

350.5 

3.1 

Y 

PMPS 

1713-1-0747 

4.57 

16.0 

28.8 

25.2 

Y 

SWHL 

1721-2457 

3.50 

47.8 

0.4 

6.8 

N 

SWIL 

1723-2837 

1.86 

19.9 

357.3 

4.2 

Y 

PMPS 

1725-38 

4.79 

158.4 

349.4 

-1.8 

Y 

PMPS 

1730-2304 

8.12 

9.6 

3.1 

6.0 

N 

SWIL 

1732-5049 

5.31 

56.8 

340.0 

-9.5 

Y 

SWIL 

1738-1-0333 

5.85 

33.8 

27.7 

17.7 

Y 

SWHL 

1741-1-1351 

3.75 

24.0 

37.9 

21.6 

Y 

SWHL 

1744-1134 

4.07 

3.1 

14.8 

9.2 

N 

SWIL 

1745-0952 

19.38 

64.5 

16.4 

9.9 

Y 

SWIL 

1748-30 

9.68 

420.2 

359.2 

-1.1 

Y 

PMPS 

1751-2857 

3.92 

42.8 

0.6 

-1.1 

Y 

PMPS 

1753-2814 

18.62 

298.4 

1.4 

-1.2 

Y 

PMPS 

1757-5322 

8.87 

30.8 

339.6 

-14.0 

Y 

SWIL 

1801-1417 

3.62 

57.2 

14.5 

4.2 

N 

PMPS 

1801-3210 

7.45 

176.7 

358.9 

-4.6 

Y 

PMPS 

1802-2124 

12.65 

149.6 

8.4 

0.6 

Y 

PMPS 

1804-2717 

9.34 

24.7 

3.5 

-2.7 

Y 

PMPS 

1813-2621 

4.43 

122.5 

5.3 

-3.9 

Y 

PMPS 

1826-24 

4.70 

81.9 

8.3 

-5.7 

Y 

PMPS 

1835-0115 

5.12 

98 

29.9 

3.0 

Y 

PMPS 

1843-1113 

1.85 

60.0 

22.1 

-3.4 

N 

PMPS 

1843-1448 

5.47 

114.6 

18.9 

-4.8 

N 

PMPS 

1853-1-1303 

4.09 

30.6 

44.9 

5.4 

Y 

PMPS 

1857-1-0943 

5.36 

13.3 

42.3 

3.1 

Y 

PMPS 

1905-1-0400 

3.78 

25.7 

38.1 

-1.3 

N 

PMPS 

1909-3744 

2.95 

10.4 

359.7 

-19.6 

Y 

SWHL 

1910-1-1256 

4.98 

38.1 

46.6 

1.8 

Y 

PMPS 

1911-1-1347 

4.63 

31.0 

47.5 

1.8 

N 

PMPS 

1918-0642 

7.65 

26.6 

30.0 

-9.1 

Y 

SWIL 

1933-6211 

3.54 

11.5 

334.4 

-28.6 

Y 

SWHL 

1934-1-1726 

4.20 

62.0 

53.2 

-1.1 

Y 

DMB 

1939-1-2134 

1.56 

71.0 

57.5 

-0.3 

N 

DMB 

2010-1323 

5.22 

22.2 

29.4 

-23.5 

Y 

SWHL 
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sar population (see , e.g., iFaucher-Giguere fc Kaspil l2006l : 
iLorimer et alJ[2006ll . The observed sample of MSPs is sub¬ 
stantially less numerous because of their generally lower lu¬ 
minosity and observational selection effects. Nevertheless, 
the large sky coverage and uniformity of the observing sys¬ 
tems of the multibeam surveys provides an excellent sample 
to begin chara cterizing their population. Recent work by 
iLorimerj (|2013lf models this population via Monte Carlo re¬ 
alizations of synthetic pulsars drawn from distribution func¬ 
tions. These synthetic populations are subsequently “ob¬ 
served” with realistic models of the surveys to produce sam¬ 
ples that can be compared with the observed data. In this 
paper, we focus on constraining the spin period distribution 
of MSPs. 

For this study, we define a MSP as a pulsar with 
P < 20 ms. Our final sample of 56 MSPs is drawn from 
detections by this survey (P MPS), the Swinburne Intermedi- 


burne High Latitude Survey (SWHL; 
the Parkes High Latitude Sur vey (PH; 


Jacobv et al. 

20071 

Bureav et al. 

2ooei 


the Perseus Arm Survey (PA; iBurgav et alJl2013r). an d the 
Deep Multibeam Survey (DMB: Ihorimer et al. 2013l i. The 
basic parameters of this sample of pulsars are summarized 
in Table [21 For the purposes of population analyses, this 
sample of pulsars is a very natural one to analyse, since the 
pulsars were found using the same telescope, receiver and 
data acquisition syst em. Since the sensitivity of this system 
is we ll understood llManchester et ^ 12001 ; iLorimer et al.l 
l2006h . our survey models are reliable. In addition, since the 
surveys were all carried out at 20 cm wavelength, we need 
not make assumptions about MSP flux density spectra in or¬ 
der to extrapolate results from surveys carried out at other 
frequencies. 


3.1 Likelihood analysis description 

Following earlier work fe.g.. ICordes &: Chernofil [l997l ~). we 
adopt a likelihood analysis to constrain the period distribu¬ 
tion. In our approach to this problem, the probability pi of 
detecting pulsar i in the sample with period Pi, and disper¬ 
sion measure DMi can be written as follows: 

Pi = f{P\a,b)V{Pi,DMi). (2) 

In this expression, as usual in statistical parlance, the 
“I” symbol denotes a conditional probability. The quantity 
/(P| a, b) represents the probability density function (PDF) 
of the period which we seek to constrain. All the models 
considered in this work can be described by two parame¬ 
ters which we refer to here generally as a and b. Specihc 
parameters will be defined below. The detectability func¬ 
tion, T>, reflects the probability of detecting the pulsar in 
one of the six surveys mentioned above and therefore appro¬ 
priately accounts for their non-uniform period sensitivity. 
Note that in this analysis, we assume that any loss of sen¬ 
sitivity due to binary motion is negligible given the signifi¬ 
cant number of ac celeration searches that have been carried 
out on these d ata ijFaulkner et al.ll20o3 : iKnispel et al.ll^lSl : 
lEatough et al.il^lSfl . 

Given a model period PDF and a detectability model 
which we describe in detail below, the likelihood function 
C{a, b) for a given combination of a and b is simply the 
product of all the 56 individual pi values. The optimal set 


of model parameters a and b are those which maximize C in 
a grid of parameter space over a and b. Once the maximum 
likelihood Tmax has been found, the marginalized PDF for 
a is then simply the distribution of T(a, b) (and vice versa 
for the PDF of b). This approach provides PDFs for a and 
b as well as a means for evaluating different period PDFs 
from the ratio of the maximum likelihoods (i.e., the Bayes 
factor, K). For example, given two model PDFs “x” and 
“y”, K — Tmax.x/Tmax,y > 1 if model X better describes 
the sample than model y. According to Jeffreys (1961), a 
Bayes factor of between 1 and 3 is deemed to be essentially 
indistinguishable from the best model while Bayes factors 
in the range 3-10 begin to favour x over y. Bayes factors 
higher than 100 decisively favour x. 


3.2 Detectability model 


The detectability of a given pulsar, U, reflects how likely it is 
to be found in the sample of 56 MSPs we will ultimately be 
applying this analysis to. Calculating 2? therefore requires 
accurately accounting for the difficulties in detecting each 
pulsar. Two approaches that can be brought to bear on this 
problem are to: (i) run large numbers of Monte Carlo simu¬ 
lations which model the detectability; (ii) develop a simple 
analytical model. After initial experimentation with the first 
approach, it became clear to us that the Monte Carlo simu¬ 
lations require a large number of assumptions and signihcant 
computational resources to carry out a sufficient number of 
realizations necessary to estimate T>. We therefore followed 
the second approach and calculate T> for each MSP in our 
sample. The essence of our approach, described in detail 
below, is to find for any line of sight the flux density distri¬ 
bution of pulsars: p{S). The detectability of a pulsar along 
this line of sight is then the fraction of such pulsars visible 
by a survey. For the pulsar, we may therefore write 


/,>(s)ds • 


(3) 


Here S'min.i is the minimum detectable survey flux density of 
this pulsar which we compute from its pulse period, disper¬ 
sion measure and pulse width (an intrinsic pulse duty cycle 
of 10% was assumed for each pulsar). The term So repre¬ 
sents the lowest flux density detectable in the survey, i.e the 
limit of S'min.i as P becomes large and DM tends to zero. 
In paper VI of this series (Lorimer et al. 2006), we gave ex¬ 
pressions for computing Smin.i and refer the reader to this 
work for details. The advantage of this approach to the prob¬ 
lem is that it does not depend on knowledge of individual 
distances to MSPs, which are uncertain. The analysis also 
does not depend on detected signal-to-noise or flux densi¬ 
ties for any individual pulsars. Instead it takes advantage of 
prior information about the pulsar population to calculate 
p{S) rigorously. As we show below, the final determination 
of the detectability function for this sample of pulsars can 
be given in terms of only two parameters which are robust 
to uncertainties inherent in the assumptions. 

The calculation of p{S) is most readily achieved from 
an application of Bayes’ theorem which implies for a given 
line of sight the following relationship between PDFs in flux 
density S and distance D-. 


p{D\S) xp{S\D)p{D). 


(4) 
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Figure 3. Detectability as a function of period (in ms) and dis¬ 
persion measure (in cm~^ pc). The data points show our esti¬ 
mates of detectability for each pulsar which we compute as de¬ 
scribed in the text. The grid shows our best approximation to 
this behaviour using the two-dimensional detectability function 
defined in equation [9] assuming the parameters ck = 10 ms and 
^ = 110 cm~^ pc. For each data point, a vertical line is drawn 
showing its distance from the best-fitting surface. 


Because the distribution we seek, p(5'), is simply p{S\D) 
marginalized over distance, we can use Eq. 0]to show that 

PiS) = J^ p{S\D)dDcxJ^ (5) 

To get expressions for p(D\S) and p(D), we use the results 
described in §3.3 of IVerbiest et al.] l|2012l '). In this work, as- 
suming a log-normal pulsar lu minosity function (see, e.g., 
iFaucher-Giguere fc Kasr)jl2006h with mean p and standard 
deviation a, it is shown that 


p{D\S) oc — exp 


1 / log S + 2 log D — p 

2 V a 


( 6 ) 


IVerbiest et ahl (l2012f ) also show that, along a line of sight 
defined by Galactic longitude I and latitude 6, an axisym- 
metric distribution of pulsars with the radial density profile 
found in paper VI leads to the result 


p{D) oc exp 

In this expression, z = Dsinb is the vertical height off the 
Galactic plane, h is the scale height of pulsars, p is a di¬ 
mensionless parameter used to scale the population over the 
Galactic disk, Rq = 8.5 kpc is the Galactocentric radius of 
the Sun and the pulsar Galactocentric radius 



o 

1 

0 ^ 

exp 

^ Ro 


(7) 


R = Rq + (D cos 6)2 — 2 J?o7Icos6cosL (8 ) 


With these analytical results it is then straightforward us¬ 
ing numerical integration of Eq. to find the appropriate 
form of p{S) for each I and b. This PDF is then numerically 
integrated according to the limit s in Eq.[3]to fi nd ’Dj . 

Following the results of iBagchi et al.l (1201 ill and 
Lorimer (2013), we adopted nominal parameter values of 
h — 500 pc, p = 5, p = —1.1 and a = 0.9. Fig. [3] shows 


the results of this calculation on our sample of 56 pulsars as 
scatter plots of D as a function of P and DM. As expected, 
T) is lower for shorter period and/or higher DM pulsars. To 
approximate this trend in the likelihood analysis, we set 

V = [1 — exp(—P/a)] exp(—DM^/2/3^), (9) 

where simple fits to the data show that a ~ 10 ms and 
/3 ~ 110 cm“® pc provide a good description of these trends, 
as shown by the smooth surface in Fig. [J] The root-mean- 
square deviation of the data from this surface is 0.1. As de¬ 
scribed in i)3.5l while the choice of duty cycle or population 
parameters assumed in the detectability analysis impacts 
the values of a and /3 somewhat, the particular values of a 
and /3 do not significantly affect our conclusions. 


3.3 Period distributions investigated 


We considered a variety of analytical functions to find the 
PDF which best describes the spin period of the MSP popu¬ 
lation. The simplest case of a uniform distribution is clearly 
not favoured by the data. In preliminary investigations we 
found Bayes factors relative to other models of order 10“^^ 
and disregarded it from further analysis. Better approxima¬ 
tions to the true period PDF can be found by considering 
functions with some well defined peak. All four functional 
forms we investigate henceforth (i.e., Gaussian, gamma, log¬ 
normal and power law distributions) require two parameters. 
In a similar way to paper VI, we refer to these parameters 
using capital letters. The Gaussian distribution has a mean 
A and standard deviation B-. 


/(P) gauss oc exp 


-jP-Af 

2B2 


( 10 ) 


The gamma distribution is parameterized by C and D: 
/(P)gamma OC exp(-P/C) (P/C)\ (11) 

The log-normal distribution is parameterized by E and P: 


/(P)lnorm OC — exp 


-(ln(P) - Ef 
2P2 


( 12 ) 


We also considered a power-law distribution parameterized 
by an exponent G and a minimum period H as follows: 


/(^)power 

= 0 forP < P, 

(13) 

/(-P)power 

oc P'^ forP > P andP < 20 ms. 

(14) 

/(^)power 

= 0forP>20ms. 

(15) 


Note that the last boundary condition simply reflects our 
definition of a MSP as a pulsar with P < 20 ms. 


3.4 Application to the observed sample 


Using the above method, we maximize C for each of the 
period distribution models. A program was writterQ to im¬ 
plement the analysis and derive marginalized PDFs of the 
resulting model parameters. For each model, we normalized 
the detection probability and period distribution such that 


P(P,DM)/(P|a,6) dP = 1. 


( 16 ) 
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Figure 4. Left and centre columns: marginalized posterior PDFs for each of the model parameters A—H (see equations 1101 — 1151 for 
definitions) obtained from the likelihood analysis described in the text. Right column: corresponding PDF for /(P) for each distribution 
when the nominal parameter values given in Table 4 are adopted. 


Table 4. Results of the likelihood analysis for the period distri¬ 
bution models considered. For each model we list the median and 
95% confidence interval on the parameters defined in equations 
1101 — 1151 along with the Bayes factor {K) computed by dividing 
the log-normal model likelihood by the likelihood of that model. 


Model 

First parameter 

Second parameter 

K 

Gaussian 

A = 

0-77o.6 

B = 5.8tj g ms 

738 

Gamma 

C = 

2.3 ±0.4 ms 

TP _ o O-I-0.4 

£/ — 

13 

Log-normal 

E = 

1.5 ±0.2 

F = 0.58l“;J2 

1 

Power-law 

G = 

-1.7 ±0.4 

H = 1.5lt° ms 

182 


This normalization ensures that the resulting likelihood val¬ 
ues can be compared with one another to compute Bayes 
factors. In the results below, we give the Bayes factors for 
the best model relative to each model under consideration. 

The results of our analysis when applied to the observed 
sample of 56 pulsars are summarized in Fig. [TJand Tabled 
Fig.|4]shows the marginalized posterior PDFs for each of the 
model parameters. Table [4] lists the 95% credible intervals 
for all the model parameters. The highest likelihood values 
were obtained for the log-normal model. The Bayes factors 
of the other models relative to this one are also given in Ta¬ 


ble |4] These results indicate that the log-normal and gamma 
distributions give by far the most plausible descriptions of 
the MSP spin period distribution. 

3.5 Testing the validity of the analysis 

Before discussing the impact of our results, it is important 
to demonstrate the reliability of the parameter estimation 
approach and its sensitivity to assumptions. To do this, we 
generated fake samples of detectable pulsars with known pe¬ 
riod distributions and passed these as input to the likelihood 
analysis. We used the psrpop softwa re packag43 introduced 
in paper VI (see also iLorimerl l2013l) to generate synthetic 
populations of MSPs for this purpose. As a starting point 
we distributed the model pulsars with model parameter val¬ 
ues of /i = 500 pc, p = 5, fj, = —1.1 and a — 0.9. For the 
period distribution, we then chose each of the four distribu¬ 
tions in turn and set the parameters A-H to be the notional 
values given in Table 0] from our analysis of the real data. 
In each simulation, we generated enough synthetic pulsars 
such that a total of 56 of them were detectable by models 

® http://psrpop.sourceforge.net 
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of the PMPS, Parkes high-latitude (PH), Perseus arm (PA), 
Deep Multibeam (DMB), Swinburne intermediate latitude 
(SWIL) and Swinburne high latitude (SWHL) surveys avail¬ 
able in psrpop. With each synthetic sample, we hrst carried 
out a detectability analysis as described in ii3.2l to deter¬ 
mine values of the detectability-model parameters a and /3 
and then applied these in our likelihood analysis. We found 
that the returned parameter values A-H from the likelihood 
analysis were entirely consistent with the input values of the 
period distribution of the parent population. In addition, 
we found that the method consistently favored the correct 
form of the input distribution by assigning it the maximum 
likelihood. For example, when we generated synthetic popu¬ 
lations assuming a log-normal distribution, we consistently 
found the Bayes factors for the log-normal likelihood model 
to be lower than the other distributions, as is seen for the ac¬ 
tual sample of MSPs. Similar results were found when other 
underlying period distributions were assumed. 

While the above results are very encouraging, they rep¬ 
resent idealized conditions in which we input the actual val¬ 
ues oi h, p, p and a into the detectability analysis to de¬ 
termine a and j3. In reality, of course, these numbers are 
not known and are only approximations to the true distri¬ 
bution of MSPs. To examine how robust the analysis is to 
changes in the assumed duty cycle, h, p, p and a, we re¬ 
peated the above procedure over a range of values to deter¬ 
mined a and /3. The ranges we explored were 5-30% duty 
cycles, 300 < h < 900 pc, 4 < p < 6, —2.5 < p < —1.5 
and 0.3 < a < 1.5. Although these led to variations in the 
detectability parameters in the ranges 2 < a < 15 ms and 
100 < P < 300 cm“® pc, we still found that the input pa¬ 
rameter distributions were recovered and that the correct 
distribution was favored. An example of this is shown in 
Fig. [5] in which we see the inferred PDFs from an analysis 
of a fake population with a log-normal period distribution. 
These results give us conhdence that our analysis on the ob¬ 
served sample of 56 MSPs is providing reliable insights into 
their underlying spin period distribution, f{P)- 


3.6 Discussion 

Based on the analysis presented in this paper, we have found 
evidence favoring the underlying spin period distribution of 
Galactic MSPs to be log-normal in form. While a gamma 
distribution is compatible with the data, it is less favoured 
than the log-normal. Uniform, power-law and Gaussian dis¬ 
tributions are decisively ruled out in our likelihood analysis 
as being good descriptions to f{P)- We note that the strong 
preference for a log-normal model found here is in contrast to 
the power-law model proposed bv ICordes fc Chernofll (Il997l ~) 
based on a much smaller sample of MSPs. While the expo¬ 
nent of our power-law model tested here (-1.7) is consistent 
with theirs, the likelihood analysis strongly favors the log¬ 
normal model. 

While our likelihood analysis weighs the different distri¬ 
butions we tested against each other, some measure of the 
absolute agreement between the log-normal model and the 
observed sample of 56 MSPs can be found by comparing 
the sample with the predicted observed period distribution 
for this model. Gombining our detectability model and log¬ 
normal period distribution, the observed period distribution 
takes the form 



Figure 6. A comparison of the sample of 56 MSPs considered 
in this paper with our best-fitting period distribution from Equa¬ 
tion 17 (solid line). 


foha{P) OC — exp 


-{\n{P)-Ef 

2^2 


1 — exp 


-P 


), ( 17 ) 


where the log-normal parameters E = 1.5 and F = 0.58 and 
the detectability parameter a = 10 ms. As can be seen from 
the comparison of this function with the binned data from 
the 56-MSP sample in Fig. the agreement is excellent, 
with the reduced value being 1.1. 

Since the sample of MSPs used in this analysis is based 
on surveys carried out a decade ago, it is useful to confront 
the distribution we obtained with the present sample of ob¬ 
jects. This is shown in cumulative form in Fig. Qwhere it is 
seen that the 95% credible region of log-normal functions we 
derive is broadly compatible with the present sample of 228 
MSPs which have been detec ted in the Parkes High Time 
Reso lution Universe Surveys llKeith et al.l lgOlOyBaxr et_al.l 
l2013l) . targeted searches of Fermi sources ( Rav et al. 20121) 
and also in s urveys at lower freq u encies with Arecib o and 
Green Bank llDeneva et al.l I2OI3I : IStovall et al.l |2014| ) . We 
note that the observed sample lies to the upper end of the 
95% credible region shown in Fig. [3 Future studies of this 
newer larger sample of MSPs should, therefore, provide more 
stringent constraints on the period distribution. 

The general agreement with our log-normal model and 
the present sample of MSPs suggests that the period- 
dependent selection effects on these “first generation” Parkes 
multibeam surveys (i.e., PMPS, PM, PA, SWIL, SWHL and 
DMB) which we model in our detectability function are 
much less severe in the present generation of MSP surveys. 


4 CONCLUSIONS 

We have presented timing models for four MSPs found as 
part of the Parkes Multibeam Pulsar survey of the Galactic 
plane. From a likelihood analysis of the sample of 56 MSPs 
detected with this earlier generation of Parkes multibeam 
pulsar surveys, we demonstrate that the underlying popu¬ 
lation of spin periods for MSPs is compatible with a log¬ 
normal distribution. When this distribution is confronted 
with more recent discoveries from other surveys, we see that 
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Figure 5. Marginalized PDFs for the log-normal parameters E and F deduced from a fake population in which the true values were 
E = 1.4 and F = 0.46 shown as dashed vertical lines. In this case, the assumed population parameters for the detectability analysis were 
intentionally biased to be /i = 900 pc, p = A, fi = —2.1 and cr = 0.5 and lead to a detectability model with a = 2.2 and /3 = 200. Even 
with such a bias, the PDFs successfully encompass the true population values and favor the log-normal model by an order of magnitude 
over the three other models. 


it is broadly consistent with the new results. It is important 
to note that the distribution we have derived here applies 
to the present-day MSP population. Although to first or¬ 
der, because of the very low spin-down rates of MSPs, the 
birth spin-perio d distribution may not be significantly dif¬ 
ferent (see, e.g., ICamilo et al.lflQM) . further investigations 
are necessary to confirm this conjecture. 

Although the true period distribution for MSPs may 
not be as simple as our analysis might initially suggest, it 
is clear that the distributions considered here all allow for 
the existence of a small fraction of pulsars with P < 1.5 ms. 
Based on the smooth curves shown in Fig. [3 the fraction of 
such pulsars in the population is around 3%. The true frac¬ 
tion could even be higher than this if we have overestimated 
the detectability of such rapidly spinning pulsars. Given our 
estimate of the analytic form of the observed period distri¬ 
bution given in Equation 17, we find that the probability of 
not detecting a pulsar in our sample of 56 MSPs with period 
P < 1.5 ms in the current sample is 99.2%. This is entirely 
consistent with the lack of such pulsars in the sample so far. 
Further discussion abou t the possibi l ity of sub-millisecond 
pulsars can be found in iLevin et al.l (l2013l i and references 
therein. 


An inspection of the current sample of MSPs shows 
no statistically significant difference between the spin pe¬ 
riod distributions of isolated objects versus binary sys¬ 
tems. A useful approach which is currently being pursued 
(P. Lazarus, private communication) is to artificially add 
short period signals to existing search data sets and directly 
test the effectiveness of pulsar search codes in recovering 
these signals. The modeling techniques presented here may 
be useful in further analyses of the MSP population which 
need to take account of the different selection biases and 
observing frequencies that have taken place since the com¬ 
pletion of the initial Parkes multibeam surveys. The tech¬ 
niques may also be applied to other population parameters 
in which selection effects may be apparent, for example the 
P — P distribution. 
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Figure 7. Cumulative distribution functions showing the ob¬ 
served sample of 56 MSPs (rightmost dashed curve), the 95% 
credible region of our best fitting log-normal model (shaded band) 
and the current sample of 206 MSPs with P < 20 ms (leftmost 
dashed curve). The deviation between the shaded band and the 
rightmost dashed curve highlights the observational selection ef¬ 
fects at short periods in our sample of 56 MSPs. 
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